## Product of the monic irreducible polynomials of degree r over F_q

from PyM import *


m = mu_moebius

def P(r,q=2):
    P = 1
    [_,X] = polynomial_ring(Z_, 'X')
    D = divisors(r)
    Dp = [d for d in D if m(d)==1]
    Dm = [d for d in D if m(d)==-1]
    for d in Dp:
        P *= (X**(q**(r//d))-X)
    for d in Dm:
        P = P/(X**(q**(r//d))-X)
    return P
    
P4 = P(4)

show(P4)

show(factor(P4))

P43 = P(4,3)

show(P43)